Role of Spatial Amplitude Fluctuations in Highly Disordered s-Wave Superconductors 
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The effect of non-magnetic impurities on 2D s-wave superconductors is studied beyond the weak disorder 
regime. Within the Bogoliubov-de Gennes (BdG) framework, the local pairing amplitude develops a broad 
distribution with significant weight near zero with increasing disorder. Surprisingly, the density of states 
continues to show a finite spectral gap. The persistence of the spectral gap at large disorder is shown to 
arise from the break up of the system into superconducting "islands". Superfiuid density and off-diagonal 
correlations show a substantial reduction at high disorder. A simple analysis of phase fluctuations about the 
highly inhomogeneous BdG state is shown to lead to a transition to a non-superconducting state. 
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The effect of strong disorder on superconductivity has 
been a subject of considerable interest, both theoreti- 
cally and experimentally for a long time. A 
generally accepted physical picture of how the supercon- 
ducting (SC) state is destroyed and the nature of the 
non-SC state has not yet emerged. Much of the theoreti- 
cal work (^pairing of exact eigenstates" or diagram- 
matics assumes that the pairing amplitude A(r) is 
uniform in space (r-independent) even for a highly dis- 
ordered SC; see however 0,^. Recent work on univer- 
sal properties at the SC-insulator transition |^ has also 
ignored amplitude fluctuations, since phase fluctuations 
are presumably responsible for critical properties. 

In this paper we consider a simple model of a 2D s- 
wave superconductor at T = in a random potential, 
defined by eqn. (1) below, and analyze it in detail within 
a Bogoliubov-deGennes (BdG) framework . Our goal 
is to see how the local pairing amplitude A(r) varies spa- 
tially in the presence of disorder, and the effect of this in- 
homogeneity on physically relevant correlation functions. 
Our results can be summarized as follows: 

(1) With increasing disorder, the distribution P(A) of the 
local paring amplitude A(r) becomes very broad, even- 
tually developing considerable weight near A « 0. 

(2) The spectral gap in the one-particle density of states 
persists even at high disorder in spite of the growing num- 
ber of sites with A(r) w 0. A detailed understanding of 
this surprising effect emerges from a study of the spatial 
variation of A(r) and of the BdG eigenfunctions. 

(3) There is substantial reduction in the superfiuid stiff- 
ness and off-diagonal correlations with increasing disor- 
der, however, the amplitude fluctuations by themselves 
cannot destroy the superconductivity. 

(4) Phase fluctuations about the inhomogeneous BdG 
state are described by a quantum XY model whose pa- 
rameters, compressibility and phase stiffness, are ob- 
tained from the BdG results. A simple analysis of this 
effective model within a self-consistent harmonic approx- 
imation leads to a transition to a non-SC state. 

We conclude with some comments on the implications 
of our results for experiments on disordered films. 



We model the 2D disordered s-wave SC by an attrac- 
tive Hubbard model with on-site disorder: 
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n = ic-\u\J2 



(1) 



^tJ2<ij>,a(4a'^jo- + h.c.) is the kinetic energy, 

(cj<j) the creation (destruction) operator for an elec- 
tron with spin cr on a site of a square lattice, t 
the near-neighbor hopping, \U\ the pairing interaction, 
"■io- = and ^ the chemical potential. The random 

potential Vi is chosen independently at each from a uni- 
form distribution [— y,F]; V thus controls the strength 
of the disorder. 

We begin by treating the spatial fluctuations of the 
pairing amplitude using the standard BdG equations pOi : 
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(2) 



where ^Unin) = -tJ2s ""(''^ + '5) + (^^ ~ A'i)^^n(j'i) and 
Au„(ri) = A(ri)u„(ri), and similarly for Vn{j^i)- Here 
S = ±x, ±y and fli /i + \U\ni/2 incorporates the site- 
dependent Hartree shift in presence of disorder. Start- 
ing with an initial guess for A(ri)'s and ili we numeri- 
cally solve for the BdG eigenvalues En and eigenvectors 
{un{ri),Vn{ri)) on a finite lattice of N sites with periodic 
boundary conditions. We then calculate the local pairing 
amplitudes and number density at T = 0, given by 



E 



(3) 



and iterate the process until self-consistency is achieved 
for Ui and A(ri) at each site, fi is determined by 
Ui = (n), where (n) is the average density. 
We have studied the model (1) for a range of parame- 
ters: 1 < \U\/t < 8 and < V/t < 12 on lattices of size 
N ~ 12 X 12 (some checks were made on 24 x 24 systems). 
We focus below on \U\/t = 2 and 4, (n) — 0.875; similar 
results are obtained for other parameters. The number of 
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iterations necessary to obtain self-consistency grows with 
disorder; we have checked that the same solution is ob- 
tained for different initial guesses. Results are averaged 
over 16-20 different realizations of the disorder. 

The distribution P(A) of local pairing amplitudes for 
\U\ = 4 is plotted in Fig. 1. For V < 0.25t, P(A) has a 
sharp peak near the V ~ BCS value of Aq — 1.36t. In 
the small V limit, pairing of exact eigenstates is justified, 
since this naturally leads to uniform A(r). However, 
this approximation fails with increasing V as P(A) be- 
comes extremely broad for V ^ t, eventually becoming 
rather skewed dX V > 2t with a large number of sites 
with A(r) « 0. 
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FIG. 1. The distribution -P(A) of the pairing amplitude 
A(ri) for different disorder strengths V . For small V , P{A) is 
peaked around the BCS Aq, but becomes increasingly broad 
at higher V, indicative of a highly inhomogeneous state. 
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FIG. 2. Density of states N{u)) for three disorder strengths 
V which show the the persistence of a spectral gap at all 
disorder; (Note the different vertical scale for each case). 

To study how the spectral gap evolves as the pairing 
amplitude becomes highly inhomogeneous, we look at the 
(disorder averaged) one-particle density of states (DOS) 
N{u;) = l/NJ2n - En), defined in terms of the BdG 



eigenvalues £"„. (Numerically, 5-functions are broadened 
into Lorentzians with a width of order spacing between 
En's). From Fig. 2 we see that with increasing disorder 
the DOS pile-up at the gap edge is progressively smeared 
out, and that states are pushed up to higher energies. 
But the most remarkable feature of Fig. 2 is the presence 
of a finite spectral gap even at high disorder. While we 
can not rule out an exponentially small tail in the low 
energy DOS from a finite system calculation, we always 
found, for each disorder realization, that the lowest BdG 
eigenvalue remains non-zero and of the order of the zero- 
disorder BCS gap; see also Fig. 3(a). We also emphasize 
that approximate treatments of the BdG equations , 
which do not treat the local amplitude fluctuations prop- 
erly, miss this remarkable feature, as do simplified models 
in which A(ri)'s are assumed to be independent random 
variables at each site. 
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FIG. 3. (a) Left panel: T = spectral gap -Egap and order 
parameter Aqp (see text) as functions of disorder V. The 
two coincide for small V but become very different at large 
disorder, (b) Right Panel: Gray-scale plot showing the spatial 
variation of A(ri) for the same disorder configuration with 
different V. Larger A(ri)'s are indicated by darker shades. 
Note the spatially correlated structures at V = 2t with "SC 
islands" separated by a "sea" of nearly vanishing A's. 

To understand the persistence of a finite spectral gap at 
high disorder, when a large fraction of the sites have near 
vanishing pairing amplitude, it is useful to study the spa- 
tial variation of of the A(ri)'s and the BdG eigenvectors 
{u„(ri)^Vn(ji)) for individual realizations of the disorder 
potential. A particularly simple picture emerges at high 
disorder: there are spatially correlated clusters of sites at 
which A(ri) is large ("SC islands"), and these are sepa- 
rated by large regions where A(ri) « (see Fig. 3(b)). 
We find that the SC islands correlate well with regions 
where the absolute magnitude of the random potential 
\Vi\ is small; deep valleys and high mountains in the po- 
tential do not allow for number fiuctuations and are thus 
not conducive to pairing. The density n{r) is also highly 
inhomogeneous, and for moderate \U\ (> At) and high 
disorder, we have found clear evidence for "particle-hole 
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mixing in real space", i.e., a spatial correlation between 
A(r) andn(r)/2[l-n(r)/2] 

At high disorder, we found that the eigenfunctions cor- 
responding to low-lying excitations live entirely on the 
SC islands (i.e., the darker regions in Fig. 3 (b)) result- 
ing in the finite spectral gap. On the other hand, regions 
where the pairing amplitude is small correspond to very 
large values of \Vi\, as explained above, and thus sup- 
port even higher energy excitations. Clearly this simple 
picture of SC islands is well defined only in the large dis- 
order regime, nevertheless, it is useful for understanding 
the spectral gap in this limit. In the opposite limit of low 
disorder, of course, the BCS-like spectral gap is obvious. 

We next turn to the question of how superconductiv- 
ity is affected in the highly inhomogeneous BdG state. 
The off-diagonal long range order parameter Aqp is 
defined by the (disorder averaged) correlation function 
(clrcliCiiCjT) ^ A^p/|C/|2 for large |r^ - rj\. From 
Fig. 3 (a) we see that Aqp is the same as the spec- 
tral gap (and both equal the uniform pairing amplitude) 
for small disorder, as expected from BCS theory. How- 
ever beyond a certain V the two quantities deviate from 
each other; in contrast to the spectral gap, the order pa- 
rameter decreases with increasing disorder; (we find that 
Aop — / (iAAP(A), i.e., the average value of the pairing 
amplitude) . 

The superfluid stiffness is given by D^/tt = 
{^kx) ~ ^xx{Qx = 0, (7j, — > 0, o; = 0). The diamagnetic 
term {—k^}, is one-half (in 2D) the kinetic energy (— /C), 
and the paramagnetic term A^x is the (disorder averaged) 
transverse current-current correlation function. We have 
also checked that the charge stiffness is equal to 
(D" is the strength of the delta-function in a{uj), and 
given in terms of A^xil = 0, a; — > 0) H). 



N=12xl2 (n) = 0.875 N=12xia (n) = 0.B75 




FIG. 4. The T = BdG superfluid stiffness D° and the 
spectral gap -Egap plotted as a function of disorder strength for 
two different values of attraction \U\. Note that irrespective of 
whether is larger than or comparable to -Egap atV = 0, the 
gap persists with increasing disorder while stiffness decreases. 



The _D?calculated within BdG theory shows a large re- 
duction |15| by two orders of magnitude with increasing 
disorder; see Fig. 4. We see that for U — — 2t, at V — 0, 
-D^ ^ ^gap, characteristic of weak coupling BCS the- 
ory, where the vanishing of the gap determines Tc, while 
for U = —4:t, D'^ and E'gap are comparable a,t V — 0, 
indicative of an intermediate coupling regime p6| where 
thermal phase fluctuations are important for determining 
Tc 0. However, for all \U\/t, we always find D° < E^^p 
at large disorder, and thus phase fluctuations have to be 
taken into account. In fact the reason why is not 
driven to zero at large V within the BdG framework is 
due to the neglect of these fluctuations. 

To make a rough estimate of the effect of phase fluc- 
tuations about the inhomogeneous BdG state we use a 
quantum XY model with an effective Hamiltonian ^ 

He = -(«/8)E,^f + P°/4)Eo-fe) cos((?, - 6^), whose 
parameters are obtained from the preceding analysis: the 
bare is the BdG superfluid stiffness and k — dn / dn \s 
the BdG compressibility. The large reduction in dn/dpi 
with disorder seen in Fig. 5 (a) can be understood qual- 
itatively at large V in terms of the charging energy of 
the SC islands. Note that, in this simplified description 
using Hg, we ignore the inhomogeneity in the local bare 
stiffness and charging energies. 

N=12xl2 (n>=0.875 U = -4t N=12xl2 (n>=0.875 U = -4l 




1 1.5 2 2.5 

V/t 

FIG. 5. (a) Left Panel: The compressibility dn/dfi as 
a function of disorder V. (b) Right Panel: The superfluid 
density Ds/vr as a function of disorder V . By including phase 
fluctuation a transition is obtained at 14 which is very close 
to that predicted by QMC 

We use a variational approximation |^8| to estimate 
the renormahzed superfluid stiffness Ds ~ Ds{k,D'^), 
by finding the best harmonic -fftriai = + 
{Ds/S)^i^-f.^^{9j — Ok)^, which describes Hg. The phase 
variables 9i are assumed to live on a lattice with lat- 
tice constant set by the BdG coherence length ^o- For 
U — —At we choose — 1-8 by demanding that 
the renormahzed Dg at = agrees with that obtained 
from quantum Monte Carlo (QMC) [|0) (Ds/n ~ 0A5t) 
for the pure case. 

We now calculate the renormahzed Ds as a function 
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of disorder, using the ^-dependent k and -0° from the 
BdG analysis as input and keeping fixed; (details will 
be presented elsewhere As shown in Fig. 5 (b), Ds 

is driven to zero beyond a critical disorder Vc which is in 
very reasonable agreement with QMC |2^] . Thus a tran- 
sition to a non-SC (insulating) state is indeed obtained 
by incorporating the effects of phase fluctuations about 
the inhomogeneous BdG state. 

We emphasize that the finite spectral gap obtained in 
the BdG analysis at large V will survive inclusion of 
phase fluctuations, since this gap is related to the in- 
homogeneous SC islands. A key question is whether the 
inhomogeneous A(r) leading to a spectral gap in the in- 
sulating state persists all the way down to \U\/t ^ 1. 
A definitive answer cannot be obtained since weak cou- 
pling BdG calculations are plagued by severe finite size 
effects It is important to note that the gap persists 
in the \U\ = 2t case (see Fig. 4(a)) which in the V = 
limit has D'^ ^ ^gapj characteristic of weak coupling 
SC. The available numerical results suggest that even for 
weak coupling, A(r) inhomogeneities are generated on 
the scale of the coherence length, which eventually show 
up as SC islands at large disorder. This would suggest 
persistence of the gap. In contrast, some tunneling ex- 
periments [Q(c) show a finite DOS N(0) with increasing 
disorder, which then points to physical effects beyond 
those in the simplest model studied here. One possibility 
is that Coulomb interactions in the presence of disorder 
lead to an effective \U\ which is r-dependent, and regions 
with \Ui\ w lead to finite A^(0). Another possibility 
is that Coulomb interactions plus disorder lead to the 
formation of local moments which are pair breaking. 

Another implication of our results for experiments is 
that SC-Insulator transitions in disordered films are often 
described in terms of two different paradigms: homoge- 
neously disordered films (driven insulating by the vanish- 
ing of the gap) and granular films (driven by vanishing 
of the phase stiffness). In our simple model, although 
the system was homogeneously disordered at the micro- 
scopic level, granular SC-like structures developed in so 
far as the pairing amplitude was concerned. It would be 
very interesting to use STM measurements to study vari- 
ations in the local density states to shed more light on 
this question. 
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